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Sufficient Conditions for the Instability of Numerical 

Integration Methods* 
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In a previous paper, a general theorem was investigated for the stability of numerical integration 
methods for the solution of systems of differential equations. In this paper, further theorems are devel- 
oped as sufficient conditions for the instability of numerical integration methods. Applying these 
theorems, the instability of known formulas are checked easily at a glance. 
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The system of m ordinary differential equations can be written in the vector form as follows: 

Y'^=F(x;Y), (1) 

with the initial value Y(a) = Fo, xe[a, b]. 
Let the vectors: 

(i) Y n be the exact solution of the following difference equation at x^Xn'. 

y„ + ,= 5; a 1 y, l+ ,.+h 2 bjr n +~ (2) 

v=-q v=-q 

(ii) Y n be the (practical) computed solution of (1) at x = x n according to: 

F n+ i= J) a v Y n +v+h J) b v F(xn+ v ;Y n+ v) + Rn+u (3) 

v=-q v=-q 

where R n +\ is the round off error, 
(iii) Y(x n ) be the exact solution of (1) at x = x n . Thus the corresponding exact formula is: 



o 



Y( Xn + h)= ^a v Y(x n + vh) + h £ b v F(x n +vh;Y(x n ^ph))-{-T k ^ (4) 



The truncation error T k+l = /x fr+ ,W +1 ^ ) 1) + 0(/i^ 2 ) , fc^l, fi k+ i * 0. 
(iv) e u be the error of the solution at the point x n , thus: 



€n=Y n -Y(Xn). (5) 
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In Karim's paper [l], 1 we discussed the following: 

1. Corresponding to the integration formula (2), there is a nonhomogeneous system of linear 
difference equations with constant coefficients for the error e n in the form: 

€ n+1 = J) d v €n+ v +hj J bvtu+v+Ri fi*0, (6) 

v=-q v=-q 

Where, in a small interval (x- q , %\): J = [j i — I is a constant matrix (i, j=l(l)m) and 

R = R n +i — Tk+\ is a nonzero constant vector, see Hamming (1959) [2]. 

2. Corresponding to the system of difference equations (6) there is a system of the character- 
istic polynomial equations: 

o 1 

P(ki)=ki- 2 a v \X-hpi ^ MT = 0, i=l{\)m (7) 

i>=-q i<=-q 

where: p/=the eigenvalues of the matrix J. 

3. The numerical integration method (2) is said to be stable according to Wilf's criterion 
(1959) [3] if the roots of (7) for alli = l(l)ra are inside the unit circle in the complex plane, 
so that the error introduced at any stage tends to decay than build up. Expanding the terms 
of (4) according to Taylor, and comparing the coefficients of h\ we get the following system 
of linear equations: 



1— j&i= 5 J(flK+^M + yiao + ^6o + y8(*+l)!M*+i, y = 0(l)*+l (8) 

yi ~ |o for ,*o< ?* = [o for 7 -4=r » = lo for ./ + *+r 



For 7 = and 7=] it follows from (8): 



£«,= !, (9) 



£i/-a„ + 2^=l. 1 10) 

v——q v——q 

REMARK: Many incorrect formulas of the form (2) are found in some references. It is necessary 
to check such cases at least by formulas (9) and (10). In many cases it is sufficient to use (9) and 
(10) for the corrections, especially in case of a sign mistake or in case of some wrong coefficients, 
otherwise it is necessary to make use of (8) fory' = 2, 3 . . ., k + 1. It is possible for example, by using 
(9) and (10) to recognize immediately that eleven of the eighteen formulas, due to H. J. Gray (1955) [4] 
(tables 4, page 14, formulas Nr. 6, 9, 10 . . ., 18) have thirteen mistakes. Using (9) and (10), the 
mistakes can be corrected. Also two formulas by W. Quade (1957) [5] (p. 160 formula (6.2') and 
p. 167 formula (9.4)) have written mistakes. 

THEOREM ll): The integration formula (2) for the solution ofm ordinary differential equations 
(1) is stable if and only if the Hermiti an forms with the coefficients 

min(r, s) _ 
A r s = ^j {Cq+i+l-r, jC q+ i +l-s. i _ C r -1, iC s -i,i}, (IX) 

1=0 



1 Figures in brackets indicate the literature references at the end of this paper. 
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for all i— 1, 2 m are positive definite 

where: 

C <r , i = -(hp i b - q H-a . q ), o- = 0(l)q+l, ai = -l, 

, f=0 trie%s . an extrapolation formula, 
b\\ if (2) is . , . r j 

[ + U an interpolation formula. 

The proof of this theorem is to be found in Karim [1|. 

THEOREM (2): //'a { , = ±1 and b_ q = — a_ q bi then the integration formula (2) is unstable. 
PROOF: Applying theorem (1), the first determinant becomes: 

Di i > = ^ffl=|C^ 1 . l | 2 -|Co. l |» 

= \hp i b 1 + a l \ 2 -\hp i b- q -r-a- q \ 2 

= \hp i b 1 -l\*-(a- q )*\hp i bi-l\* = 0. 

Thus, it is sufficient that, the integration formula (2) is unstable. 

l (q-s)/2 [ Q even 

Theorem (3): //^ (— l)"b„ = and^ a_ q+2 p=s = j . ifqis ... 

then the integration formula (2) is unstable. 
Proof: 

1. If (j is an even number, we get 

2 r/ 9+21/ H" J)" "+ 2 " ' [+1 X" 9 +2t; ~ 2" 9+2^-1 [ = 0. 

Applying (9) it follows: 

£ (-l)'a„=-l. 

l'=-f/ 

Using the first assumption of the theorem we get: 

2 (-1 )'*«„+% J) (-1)"6,. = -1 (12) 

p=-g p=-q 

2. If g is an odd number, we get: 

r(«-D/2 (9-O/2 -I r(«-D/2 (g+D/2 <] 

Also 

(g+0/2 (r/-l)/2 

^j Cl-q+2v-l ~~~ 2j a- q+ 2u~—l, 
v=\ v=Q 



this means V (— l) v a v = — 1. 

Applying the first assumption of the theorem we get: 



2 {-l)'a.+ hpi X (-l)'6r— 1. (12') 
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According to (12) and (12'), Xf = — 1 is a characteristic root of the equation (7), whether q is an odd 

or even number. Thus, due to Wilf's definition, the integration formula (2) is unstable. 
i 
Theorem (4): // j b„=0, (13) 

then: 

(i) The system of the nonhomogeneous difference equations (6) has a nonconstant particular 

solution, 
(ii) The integration method (2) is unstable. 
Proof: 
(i) To find a particular solution of (6), let 

e v = e= a nonzero constant vector, v = — q(l)l. (14) 

Thus we get from (6): 

(l- J) a v )e = hje 2 6 "+*- 

\ v~-Q I v--q 

R = 0. (15) 



Using (9) and (13) it follows: 



This means that the assumption (14) contradicts (15) with reference to (6). Put e v =ve in (6), where 
e is a nonzero constant vector. Using (9) and (13) we get: 

(l- 2 *, 'A * = h J 2 h *- v + R- 

\ i>=-q I v=—q 

Hence assuming that 2 b v * v 4 s 0, det J =¥ 0, and applying (10) and (13), it can be proved that: 

v=-q 

e = -(h J b,-^ j-tR. 

\ v=-q f 

Thus the equations system (6) has a particular solution which is equal to 

€n = -n(h 2 bu-v) J' l R. 

1 

In the case of 2 ^ ' v = ®i then e M is a polynomial in n of degree greater than 1. 

v=-q 

(ii) Using (9) and (13) we have: 

o 1 

1— 2 a-v — hpi 2 b v = Q. 

i'——q v=—q 

This means X,= l is a characteristic root of the eq (7). Thus, and in addition of (i), the stability 
condition is not satisfied and the integration formula (2) is unstable. 

Illustrative Examples 

The instability of the following integration formulas may be easily checked at a glance: 

(1) Y n+1 = Y n - l + 2hY^ T*=\h?Y"(xn) + 0{h A ) 

122 



(2) y. + i— y,_i + k„_, + f„ +&(- f;,_, + y;,) , t, = h*Y" (*„) + o(h<) 

0) F„ + ,=-y„- 2 + r„-i + r„+2M- r»_,+ n) , r 4 =i w^*.) + o(A 5 ) 

(4) y H+ ,=-r„_,+2K„+| (-F;,_, + y;, +1 ), r^-^^i™ (*„)+<)( A 5 ) 

4A 14 

(5) F„ +1 = y„- 3 + T (2FU-f;,_,+2f;,), n^ iwm +o(A 6 ) 

(6) f„ +1 = f„_ 3 +f h _ 2 -f„+3/ 1 (f;,_ 2 + f;,), n=^^i A ^„) + o(/ l < i ) 



(7) F„ + , = F„_ 2 + 9F„_, - 9F„ + 6A( F,',_, + F/,) , r, = jo * 5yV ^ + 0( - h ^ 

(8) f„ +i =-f„_ 2 + f h - 1 + f,+| (-f;,_ 2 -3f;,__, + 3f;, + f;, + ,), r«=-^/i 6 y v, K)+o(A 7 ) 

? 41 

(9) f„ j , = f„_,+ T q/ 1 (iif;^-i4f;,_ :i +26f;,„ 2 T 7 =j^h j r /u (x„)+o(h*) 

-uyi^ + uy;,). 

(10) F„ fl =2^8 (-^8F,-., + 243F„_ 4 + 250F„_ 2 + 243F„) 7k=_ 4^ ^^ m (^) + 0(A 9 ). 

+^/>(-8F;,_ s -27F;,_ 4 +27r;,+8F;, +r ), 



According to theorem (2), the formulas (1), (5), (6), (7), and (9) are unstable. According to theorem 
(4), the formulas (2), (3), (4), (8), and (10) are unstable. Moreover, according to theorem (3), the 
formula (2) is unstable. Thus all the given formulas (1), (2), . . ., (10) are unstable. 

Note: The formulas (1), (4), (6), (8), and (10) are given by Quade (1957) [5], (5) is given by 
Milne (1949) [61, (3) is given by Hamming (1959) [2], (9) is given by Wilf (1960) [7], (2) and (7) are 
given by the author. 
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